/*! 
 ***************************************************************************
 * \file   convection2.cpp
 * Forced Convection-Diffusion in 2D.
 *
 ***************************************************************************
 * \author    Luis M. de la Cruz  [ Sat May 23 12:06:36 CDT 2009 ]
 ***************************************************************************
 * \par Description
 * This code solves the next equation:
 * \f[ 
 * \frac{\partial T}{\partial t} + 
 * u \frac{\partial T}{\partial x} + v \frac{\partial T}{\partial y}  = 
 * \frac{\partial T}{\partial x^2} + \frac{\partial T}{\partial x^2}
 * \f]. 
 * \par
 * This equation is solved in a unit square \f$ x, y \in [0,1] \f$. 
 * The velocity is prescribed as follows:
 * \f[
 * u = -A \cos(\pi y) \sin(\pi \lambda x / L_x );
 * \f]
 * \f[
 * v = \frac{A \lambda}{L_x} \sin(\pi y) \cos(\pi \lambda x / L_x);
 * \f]
 * The boundary conditions are as shown in the next figure:
 *
 * \image html convection.png "Unit square in 2D" width=10cm
 * \image latex convection.eps "Unit square in 2D" width=10cm
 *
 * The equation can be solved using Upwind, Central Average (CDS) or QUICK 
 * schemes for the approximation of convective terms. CDS and QUICK can
 * be used with \c CDS_CoDi and \c QUICK_CoDi, or with \c CDS_Hay and 
 * \c QUICK_Hay. The latter two are the deferred implementation from
 * Hayase et al. 1992, which are more stable.
 *
 * The results are as presented in the next figure:
 *
 * \image html  convForced2D01.png "Final result" width=10cm 
 * \image latex  convForced2D01.eps "Final result" width=10cm 
 *
 *  \par Compiling and running
 *  Modify the variables BASE and BLITZ in the file \c tuna-cfd-rules.in 
 *  according to your installation and then type the next commands:
 * \par
   \verbatim
   % make
   % ./convection2 \endverbatim
 * 
 ***************************************************************************
 */

#include <iostream>
#include <cmath>
#include "Meshes/Uniform.hpp"
#include "Equations/ScalarEquation.hpp"

#ifdef WITH_CDS_CoDi
#include "Schemes/CDS_CoDi.hpp"
#endif

#ifdef WITH_CDS_Hay
#include "Schemes/CDS_Hay.hpp"
#endif

#ifdef WITH_Upwind_CoDi
#include "Schemes/Upwind_CoDi.hpp"
#endif

#ifdef WITH_Quick_CoDi
#include "Schemes/Quick_CoDi.hpp"
#endif

#ifdef WITH_Quick_Hay
#include "Schemes/Quick_Hay.hpp"
#endif

#include "Solvers/TDMA.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"
using namespace Tuna;
  
// Scalar field in one dimension
typedef TunaArray<double,2>::huge ScalarField2D;
typedef double prec_t;

template<class sfield>
void initialVelocity(sfield &u, sfield &v, prec_t dx, prec_t dy,
		     prec_t A, prec_t l_x, prec_t l)
{
    prec_t A2 = A * l / l_x;
    const int bi = u.lbound(firstDim) + 1, 
	ei = u.ubound(firstDim) - 1,
	bj = u.lbound(secondDim) + 1, 
	ej = u.ubound(secondDim) - 1;
    
    for(int i = bi; i <= ei; ++i)
	for(int j = bj; j <= ej; ++j) {
	    u(i,j) = -A * cos (PI * dy * (j-0.0)) * 
		sin (PI * l * dx * (i-0.0) / l_x );    
	    v(i,j) = A2 * sin (PI * dy * (j-0.0)) * 
		cos (PI * l * dx * (i-0.0) / l_x);
	}
}

int main() 
{

#ifdef WITH_GNUPLOT
    GNUplot plotter;
    int stop;
    plotter("set grid");
    plotter("set contour");
    plotter("set cntrparam level discrete 0.4, 0.3, 0.2, 0.1, 0, -0.1, -0.2, -0.3, -0.4");
    plotter("set view 0,0");
    plotter("set size 0.73,1");
    plotter("unset surface");
#endif

    prec_t length_x, length_y;
    int num_nodes_x, num_nodes_y;
    prec_t right_wall = -0.5, left_wall = 0.5;
    prec_t tolerance, dt, error = 1.0, residual;
    prec_t lambda, amplitude;
    int max_iter, nmax, frequency, iteration = 0;
    timer time;

// Reading data from the input file
    std::ifstream input_cin ("input2");
    input_cin >> length_x
	      >> length_y
	      >> num_nodes_x 
	      >> num_nodes_y 
	      >> dt
	      >> nmax
	      >> tolerance
	      >> max_iter
	      >> lambda
	      >> amplitude
	      >> frequency;
    input_cin.close();

    StructuredMesh<Uniform<prec_t, 2> > mesh(length_x, num_nodes_x,
					     length_y, num_nodes_y);
    mesh.print();
    prec_t dx = mesh.getDelta(X);
    prec_t dy = mesh.getDelta(Y);
    int num_vols_x = mesh.getExtentVolumes(X);
    int num_vols_y = mesh.getExtentVolumes(Y);

    ScalarField2D T(mesh.getExtentVolumes());
    ScalarField2D u(mesh.getExtentVolumes());
    ScalarField2D v(mesh.getExtentVolumes());
    ScalarField2D us(num_nodes_x, num_vols_y); // staggered u-velocity 
    ScalarField2D vs(num_vols_x, num_nodes_y); // staggered v-velocity

    ScalarField2D Tn(mesh.getExtentNodes());
    ScalarField2D un(mesh.getExtentNodes());
    ScalarField2D vn(mesh.getExtentNodes());

    Range all = Range::all();
    T(T.lbound(firstDim), all) = left_wall;   
    T(T.ubound(firstDim), all) = right_wall;

    T(all, T.lbound(secondDim)) = left_wall;   
    T(all, T.ubound(secondDim)) = right_wall;   

    initialVelocity(un, vn, dx, dy, amplitude, length_x, lambda);
    NumUtils::staggerFromNodesX(us, un);
    NumUtils::staggerFromNodesY(vs, vn);

#ifdef WITH_DX
    InOut::writeToFile_DX(un, vn, 0, "./Data2/velc.", dx, dy);
#endif

    SparseMatrix< Diagonal<prec_t, 2 > > A(num_nodes_x, num_nodes_y); 

    ScalarField2D              b(num_nodes_x, num_nodes_y);

#ifdef WITH_CDS_CoDi
    ScalarEquation<CDS_CoDi<double,2> > conv_diff(T, A, b, mesh.getDeltas()); 
#endif

#ifdef WITH_CDS_Hay
    ScalarEquation<CDS_Hay<double,2> > conv_diff(T, A, b, mesh.getDeltas()); 
#endif

#ifdef WITH_Upwind_CoDi
    ScalarEquation<Upwind_CoDi<double,2> > conv_diff(T, A, b, mesh.getDeltas());
#endif

#ifdef WITH_Quick_CoDi
    ScalarEquation<Quick_CoDi<double,2> > conv_diff(T, A, b, mesh.getDeltas());    
#endif

#ifdef WITH_Quick_Hay
    ScalarEquation<Quick_Hay<double,2> > conv_diff(T, A, b, mesh.getDeltas());    
#endif

    conv_diff.setDeltaTime(dt);
    conv_diff.setGamma(1.0);
    conv_diff.setNeumann(TOP_WALL);
    conv_diff.setNeumann(BOTTOM_WALL);
    //    conv_diff.setDirichlet(TOP_WALL);
    //    conv_diff.setDirichlet(BOTTOM_WALL);
    conv_diff.setDirichlet(LEFT_WALL);
    conv_diff.setDirichlet(RIGHT_WALL);
    conv_diff.setUvelocity(us);
    conv_diff.setVvelocity(vs);
    conv_diff.print();

    ofstream error_fp("./Data2/error"), residual_fp("./Data2/residual"), 
      nusselt_fp("./Data2/nusselt");

    double total_time = 0.0;
    while ( (error > tolerance) && (iteration <= nmax) ) {

        time.tic();
        conv_diff.calcCoefficients();
	Solver::TDMA2DY(conv_diff, tolerance, max_iter);
	residual = conv_diff.calcResidual();
	error = conv_diff.calcErrorL1();
	conv_diff.update();
	total_time += time.toc();

        nusselt_fp << iteration  << "\t"
		   << NumUtils::calcNusseltVert(T,dx*0.5,dy,0) << endl;
	
	error_fp << iteration  << "\t" << error << endl;
	residual_fp << iteration  << "\t" << residual << endl;

	cout << "\n ---- Iteration = " << ++iteration
	     << " Residual = " << residual
	     << " Error = " << error
	     << " Total time = " << total_time
	     << endl;

	if (!(iteration % frequency)) {
	    NumUtils::interpolateToNodes(Tn, T);

#ifdef WITH_DX
	  InOut::writeToFile_DX(Tn, iteration, "./Data2/temp.", dx, dy);
#endif

#ifdef WITH_GNUPLOT
	    InOut::writeToFile(Tn, 1, "./Data2/temp.", dx, dy);
	    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
	    std::cin >> stop;	
	    plotter("splot \"Data2/temp.00001\" w l");
#endif

	}	
    }

#ifdef WITH_DX
    std::cout << "\n\n NOTE: The data were written in the [./Data2] directory in DX format";
#endif

    return 0;
}







